Kinetic Monte Carlo simulation of faceted islands in heteroepitaxy using multi-state 

lattice model 
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A solid-on-solid model is generalized to study the formation of Ge pyramid islands bounded by 
(105) facets on Si(lOO) substrates in two dimensions. Each atomic column is not only characterized 
by the local surface height but also by two deformation state variables dictating the local surface 
tilt and vertical extension. These deformations phenomenologically model surface reconstructions 
in (105) facets and enable the formation of islands which better resemble faceted pyramids. We 
demonstrate the model by application to a kinetic limited growth regime. We observe significantly 
reduced growth rates after faceting and a continuous nucleation of new islands until overcrowding 
occurs. 

PACS numbers: 68.65.Hb, 81.16.Dn, 81.16.Rf 



I. INTRODUCTION 

Strain induced self-assembly of three dimensional (3D) 
islands in heteroepitaxy have been attracting much re- 
search interest because of the rich physics involved and 
their potential applications as quantum dots in optoelec- 
tronic devices [1113, S|- A widely studied system is Ge 
deposited on Si(lOO) substrates with a 4% lattice misfit. 
Relatively flat islands in the form of stepped mounds with 
unfaceted sidewalls called pre-pyramids start to emerge 
at 3 monolayers (MLs) of Ge coverage [3, Further 
deposition leads to pyramids or rectangular-based huts 
bounded by (105) facet planes. Deposition temperatures 
lower than 500°C generally favors rectangular huts [1, 0| 
while higher temperature often leads to pyramids [1] . Af- 
ter still further deposition or annealing, pyramids can 
grow into dome islands bounded mainly by steeper (113) 
facets [1,113. 

(105) facets on pyramids and huts have been found to 
be extraordinarily stable and atomically flat from flrst 
principle calculations [lllj [H) iHj At low temper- 

ature, surface steps on (105) facets are rarely observed 
[3 . They are however present at higher temperature and 
the bunching of them are observed to be important to 
the morphological evolution [T5| . The structures, ener- 
gies and dynamics of these steps have been studied using 
first-principles calculations [ly|. Also, the edge energies 
of a (105) faceted ridge have been estimated using molec- 
ular dynamics simulations based on empirical potentials 

[ill]- 

Large scale simulations of the formation of 3D islands 
is possible using kinetic Monte Carlo (KMC) methods 
based on lattice models [11 [H, S HI [H, M, Hi [H, 
m, ll^l ■ The simulations are computationally very inten- 
sive due to the long-range nature of elastic interactions. 
Elastic forces can be accounted for accurately and effi- 
ciently using advanced alg orithms so that simulations in 
2D [llill and 3D [Hill [H, Hi] with respectively large 
and moderate system sizes are possible. Using more ap- 
proximate forms of the elastic interactions, larger systems 
in 3D can also be studied [20, ^3 ■ 



KMC studies on strained layers are generally based on 
square or cubic lattices for simplicity. Strain induced is- 
lands or pits are readily generated but their sidewalls are 
almost vertical Il8l. l24'. 261 or at an of inclination of about 
45 degrees [Til. l2ll |23, 27j depending on the details of 
the bond energies or additional constraints used. These 
inclinations are much steeper than 11° and 26° for the 
(105) and (113) facets respectively. The realistic facets 
however are of rather low-symmetry and in general are 
not favored energetically in lattice models. The discrep- 
ancy results in strain distributions considerably different 
from the realistic ones and may probably lead to qualita- 
tively different growth modes in certain situations. Fur- 
thermore, the surface energy of the island sidewalls from 
existing KMC models are not independently adjustable 
and there is no simple approach to incorporate for ex- 
ample the extraordinary stability of certain facets. In 
addition, with only one favored sidewall slope in a given 
model, only one type of island can be simulated so that 
studying the pyramid to dome transition for instance is 
impossible. 

In this work, we extend the convectional ball and 
spring lattice model for KMC simulation of heteroepi- 
taxial solids in 2D by allowing specific geometrical defor- 
mation states of the surface atoms. These deformations 
phenomenologically represent surface reconstructions on 
(105) facets. We show computationally that this new 
multi-state model leads to the formation of faceted is- 
lands. Examples of qualitative differences in the growth 
dynamics between faceted and unfaceted islands are ex- 
plained. 



II. BALL AND SPRING LATTICE MODEL 

We first explain the conventional square lattice model 
of clastic solids in 2D while further extensions will be in- 
troduced in the next section. Every atom is associated 
with a lattice sites and are connected to nearest and next 
nearest neighbors by elastic springs. Solid-on-solid con- 
ditions are assumed. We follow the model parameters 



2 



used in Ref. [2l| unless otherwise stated to approximate 
the widely studied Ge/Si(001) system. We assume a sub- 
strate lattice constant = 2. 71 5 A so that gives the 
correct atomic volume in crystalline silicon. The lattice 
misfit e — {af~as)/af equals 4% where a/ is the lattice 
constant of the film. Nearest and next nearest neigh- 
boring atoms are directly connected by elastic springs 
with force constants fc^v = 13.85ey/a^ and kj^N = /2 
respectively. The elastic couplings of adatoms with the 
rest of the system are weak and are completely neglected 
for better computational efhciency. In this model, surface 
steps have a particularly high tendency to bunch together 
under strain presumably due to the much weaker entropic 
surface step repulsion in 2D. We hence forbid double sur- 
face steps as well as adjacent single surface steps of the 
same direction so that the steepest surface slope allowed 
is 1/2. 

The KMC approach simulates the morphological evo- 
lution by explicitly considering the diffusion of surface 
atoms. Every topmost atom m on the film can hop to 
a nearby site with a hopping rate T(rn) following an Ar- 
rhenius form: 



r(m) — Rq exp 



knT 



(1) 



where Um is the number of nearest and next nearest 
neighbors of atom m. We have assumed an identical 
nearest and next nearest neighbor bond strength 7. We 
put 7 = 0.5eV, sHghtly larger than the value in Ref. 
[2l| so that the energy costs of stepped mounds become 
slightly higher. The energy AEs{m) is the difference in 
the strain energy Eg of the whole lattice at mechanical 
equilibrium with or without the atom m. Due to the 
long-range nature of elastic interactions, its efficient cal- 
culation is highly nontrivial and we handle it using a 
Green's function method together with a super-particle 
approach explained in Refs. [U, [2^ [2^. In addition, 
Eq = 87 — 0.67eV, where 0.67eV is the adatom diffu- 
sion barrier on the (100) plane. To speed up the simula- 
tions, long jumps are allowed so that a hopping atom will 
jump directly to another random topmost site at most 
Smax = 8 columns away with equal probability. Then, 
Rq = 2Do/{asasf with Do = 3.83 x lO^^^^g-i ^j^j 
= \{s„iax + l)(2s,„Q3; + 1). This gives the appropriate 
adatom diffusion coefficient for silicon (100). 



III. MULTI-STATE LATTICE MODEL WITH 
SURFACE DEFORMATION 



although correlations between misfit strain and surface 
reconstruction are known to exist, [l^, [13, [13] ■ In the 
following calculation of the local deformation energies, 
we hence neglect lattice misfit and express all lengths in 
unit of lattice constant. The subsequent calculation of 
the misfit strain energy term is identical to that outlined 
in Sec. El 



(a) 



■■■li 



FIG. 1: A faceted island from a small scale simulation using 
the multi-state model (a) and a magnification of part of the 
surface containing a (105) surface step between the third and 
the fourth columns (b). Deformed film atoms, undeformed 
film atoms and all substrate atoms are shaded in red, light 
blue and dark blue respectively. In (b), the tilt variable Ui 
is i for all columns, while the extension variable Ki from left 



to right equals 0, i |, -i, 0, 



0, i, 1,-1, -i and 



We first show an example of a faceted island from a 
small scale simulation in Fig. [TJa). Figure[ljb) magnifies 
part of the surface. It shows how the surface deforma- 
tion smooths out the (100) steps of the original stepped 
mound and turn the sidewalls into atomically flat effec- 
tive (105) facets with slopes ±1/5. An example of a sur- 
face step on the (105) surface is also shown and will be 
explained later. In the absence of deformation, an atom 
is represented by a unit square. An integer hi denotes 
the surface height at column i. We assume that a top- 
most atom in the film surface or in an exposed region of 
the substrate can be deformed into a trapezoid charac- 
terized by two new deformation state variables, namely 
a tilt variable ui and an extension variable Ki . We put 



To effectively model (105) facets, which are more pre- 
cisely (15) surfaces in 2D, we introduce additional de- 
grees of freedom representing local deformations to all 
topmost atoms. They phenomenologically accounts for 
the surface rebonding or reconstruction states on a (105) 
faceted region [ll|. For efficient computation, these de- 
formations localized to individual surface atoms are as- 
sumed to be completely independent of the lattice misfit. 



1 

CTi = 0, -, 
5 



or 



(2) 



which gives the slope of the upper surface of the deformed 
atom. The values ai = ±1/5 enable the formation of 
the (105) facets in both directions. As shown in Fig. 
[IJb), attaining a flat (105) faceted region further requires 
properly coordinated vertical stretching or compression 
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of the topmost atom by Ki which is given by 



,0 i, or 



for ai — 
for CTi = ±T 



(3) 



The ith atomic column hence can be rectangular or trape- 
zoidal with the left and right edges of heights hf and 
given by 



ht = h, 
h'' = h,. 



2 

<Jj_ 
2 



(4) 
(5) 



A surface step in between the «th and the (i+l)th column 
has a step height 5i defined as 



(6) 



For simplicity, we have measured step heights as pro- 
jected along the lattice axis rather than the surface nor- 
mals. Note that single steps on (100) and (105) surfaces 
have very different heights of 1 and 1/5 respectively in 
our model. 

We will next explain the energy cost of the local defor- 
mation of the surface atoms. Values of the energy param- 
eters to be introduced are chosen phenonienologicaUy to 
provide morphologies best compared with observations. 
Similar to the original lattice model |2l|, although we 
believe that our parameters are within physically accept- 
able ranges, this model being in 2D is not realistic enough 
to appl y d irectly parameters from first principle studies 
[r^.ll3.Tl^ in general. Furthermore, we have found from 
numerous exploratory simulations that only a rather lim- 
ited and specific range of parameters provides reasonable 
morphologies under a wide range of relevant growth con- 
ditions. The constraints on our parameters hence may 
also shed light on how the morphologies reveals certain 
features on the microscopic details of the surface and this 
will be discussed further. 

The hopping rate of a topmost atom m in Eq. ([T]) is 
generalized to 



r(m) = i?o Gxp 



AEb{m) + AEs{m) 



(7) 



where 



-7 — 0.67 eV. The misfit strain energy term 
AE's (m) is defined similarly as before and its calculation 
is assumed to be completely independent of the local sur- 
face deformation. The surface energy term A£'f,(m) de- 
notes the change in the bond energy of the whole sur- 
face when the site is occupied versus unoccupied. More 
precisely, surface energy is defined relative to that of a 
flat (100) surface as 

Eb = [vi'^i) + ^{f^h o-j+i) -I- ujiSi, (Ti, (Ti+i)] (8) 



Here, ri(±l/5) = 5 meV is the formation energy per site 
of the (105) facet and 77(0) = for the (100) region. Also, 



i/(cri, cTi+i) = 0.35 eV denotes the interface energy at the 
boundary of a facet where ai ^ (Ji+i and it is zero other- 
wise. It dictates the energy barrier of facet nucleation. If 
we choose a larger value of 77(±l/5), the (105) facet can 
become unstable. A negative value of ri{±l/5) has been 
suggested [l3| corresponding to extremely stable (105) 
facets. However, this is not acceptable as island sizes 
from such simulations are then dominated by v closely 
related to the edge energy in Ref. [l3| but is practically 
independent of the lattice misfit. 

The last term in Eq. ([5]) represents the energy of a 
surface step. On a (100) region with ai = ai+i = 0, it is 
defined as 



W((5i, (Ti, (Ti+l) = 



(9) 



where the step height 5i defined in Eq. ([S]) is a integer. 
This results from simple bond counting noting that two 
single steps are created by breaking one nearest neigh- 
boring bond of strength 7. Noting also that a bulk atom 
has a bond energy of —47, Eq. reduces exactly to 

Eq. ([1]) so that the (100) regions in the multi-state model 
behaves identically to the basic model in Sec. |TT1 Outside 
of a (100) region (i.e. ai or (Ti+i ^ 0) we put 



uj{5i 



a„ a,+i) = /3io5 (1 + X - xe^'"^^') + 



(10) 

for 5i > 1/5 and it is zero otherwise. This expres- 
sion gives an energy /3io5 for a single step with height 
Si = 1/5 on a (105) region. It is known that incomplete 
(105) facets can be practically absent at low temperature 
around 450°C 0] but are observable at 550°C [3. We 
reproduce this feature in our model by taking a relatively 
large value of /3io5 = 0.3 eV. From Eq. (fTU|) . the step en- 
ergy per unit height of a multiple step approaches 7n/2 
identical to that for a step on a (100) facet. This also re- 
duces the energy of an adatom on a (105) surface which is 
bounded by two unit steps to a more acceptable but still 
very large value of 1.3 eV. The parameter x determines 
the energy of multiple steps of intermediate heights. We 
put X = 0.5 allowing a slight tendency of step bunching 

In KMC simulation using this multi-state model, the 
atomic hopping events are randomly sampled and simu- 
lated according to the rates T{m) in Eq. ([7]). We assume 
that the deformation state variables ai and Hi at every 
column are unchanged after an atomic hop, i.e. the defor- 
mation state is attached to the column rather than to the 
hopping atom. Deposition of an atom also increases the 
column height by unity without altering the deformation 
state. After every period t, the deformation state for a 
set of columns will be updated. Specifically, to facilitate 
program parallelization, we adopt a sublattice updating 
scheme in which the deformation states at all odd (even) 
lattice sites will be updated at every odd (even) updating 
event. When column i is to be updated, the variables ai 
and Ki are re-sampled from the allowed set of 11 possible 
combinations using a heat bath algorithm based on the 
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relative probability exp{—Eb/kT). We take r = 2 /Tad 
where Tad is the adatom hopping rate on a (100) sur- 
face easily calculable from Eq. Ilj. This is the highest 
possible rate without increasing significantly the overall 
execution time of our program. Local changes in the sur- 
face reconstruction states are most likely a fast process 
compared with atomic hopping. We have checked that 
our deformation state updating rate is indeed sufhciently 
fast so that decreasing r gives no observable difference 
to our results. Our model follows detailed balance which 
allows us to confirm the reliability of our software imple- 
mentation using a Boltzmann's distribution test [25| . 



IV. RESULTS 

Using both the conventional ball and spring lattice 
model and the multi-state lattice model with surface de- 
formation explained in Sees. |TT] and IIIIl we have sim- 
ulated the self-assembly of strained islands in 2D. A 
substrate of size 1024 x 1024 (width x depth) is used. 
We take a temperature 450° and a deposition rate 0.1 
ML/s. The conventional and the multi-state models lead 
to islands with unfaceted and faceted sidewalls respec- 
tively. For convenience, we refer to them as unfaceted 
and faceted islands. 

Figure El^a) shows the evolution of unfaceted islands 
from a typical run using the conventional model during 
deposition of up to 6 MLs of film material on to an ini- 
tially flat substrate. Unstable shallow stepped mounds 
develop at very early stage. After depositing about 2 
MLs, some stepped mounds have attained steeper side- 
walls and become more stable. At about 4 MLs, they 
have generally attained the steepest possible slope of 1/2 
allowed in our model. As observed in this and other 
similar runs, there is a rather well defined island nucle- 
ation period and no new island emerges after some larger 
islands are established. We also observe that some rel- 
atively mature islands eventually decay and vanish, in- 
dicating a ripening process. The existence of a finite 
nucleation period followed by ripening is consistent with 
previous KMC simulations [27] as well as continuum sim- 
ulations [29|, (3^. It may also have some experimental 
relevance at higher temperature although the pyramid 
to dome transition and alloying between the film and 
substrate atoms [To"] add further complications. 

Analogous evolution of faceted islands simulated us- 
ing the multi-state model with surface deformation is 
shown in Fig. (Hb). Small highly unstable (105) faceted 
regions with deformed surface atoms begin to appear at 
a coverage of about 0.5 MLs. Relatively stable (105) 
faceted islands emerges at about 1 ML. These islands 
develops from the larger ones of the stepped mounds. 
Faceted regions nucleate on either side of the mounds in- 
dependently so that half faceted asymmetric islands exist 
during the course of development. Islands also often go 
through an truncated pyramid stage 5] with unfaceted 
tops before finally becoming fully developed pyramids. 
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FIG. 2: Snapshots of surfaces showing the development of (a) 
unfaceted and (b) faceted islands simulated respectively us- 
ing the conventional model and the multi-state model. (105) 
faceted regions are shaded in red. Each successive profile is 
displaced by -1-5 vertically for clarity and corresponds to the 
deposition of a further 1/4 MLs up to a total of 6 MLs. 



Some faceted islands may occasionally decay partially 
or even completely back to unfaceted stepped mounds, 
but the larger ones are much more stable. On the other 
hand, some stepped mounds may happen to get faceted 
at rather small sizes while slightly larger ones can remain 
unfaceted for long periods. Therefore, the faceting pro- 
cess in our current model is strongly affected by both the 
energetics and the kinetics. 

At this low growth temperature of 450°, surface steps 
on a (105) facet is rare as explained in Sec. IIIIl Further 
growth of faceted islands by step flow is hence kinetic lim- 
ited [1, 0| ■ It can be observed from Fig. Hfb) that island 
growth rates drop dramatically once becoming faceted. 
Their sizes occasionally jump up rapidly only when parts 
of the sidewalls become temporarily unfaceted due to 
thermal excitation. Since developed islands are poor ab- 
sorber of newly deposited atoms, new islands continue 
to nucleate until the substrate is crowded with islands. 
Kinetic limited growth and continuous island nucleation 
have not been reported previously in KMC or continuum 
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simulations in our knowledge. More importantly, depo- 
sition experiments at 550° C do indicate slower growth 
of matured islands and a continuous island nucleation 
growth mode [8]. Deposition at lower temperature how- 
ever leads to huts [y, 0| which may share some related 
characteristics but are more complicated. 
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FIG. 3: Plot of island size against nominal film thickness h 
for unfaceted (a) and faceted (b) islands 

For more quantitative analysis, we define an island as 
one in which each of the constituent columns must be at 
least 4 atoms tall. All islands can then be automatically 
identified. Figure [3] traces the size evolution against the 
nominal film thickness h of every island in Fig. [5] once 
they have attained a size of at least 150 atoms. Islands 
from another similar run are also included in Fig.[3Ub) to 
provide additional examples. From Fig. [SJa), unfaceted 
islands beyond a certain size in general grow steadily with 
its own characteristic rates which are expected to depend 
mainly on the sizes of their adatom capture zones. Small 
islands decay and vanish. In contrast, from Fig. EJb), 
there is in general an initial period of rapid island growth 
followed by much slower growth after faceting. Once 
faceted, their sizes remain nearly constant except at oc- 
casional jumps associated with temporary partial decay 
of the facets as described above. 

To obtain more statistics, we have repeated each sim- 
ulation 200 times. Figure [3] plots the average number of 
islands of size 150 or larger on the 1024 atoms wide sub- 
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FIG. 4: Plot of the average number of islands on a substrate 
of 1024 atoms wide against nominal film thickness h 



strate used. Smaller islands are excluded because they 
are highly unstable. For unfaceted islands, their number 
first increases indicating a period of active nucleation at 
coverage from about 1 to 2.5 MLs. It then declines but 
at a very slowly rate indicating rather inefficient coarsen- 
ing during growth. In contrast, faceted islands steadily 
increase in number for coverage up to about 5 MLs due 
to continuous nucleation. Beyond 5 MLs, the substrate 
is crowded with islands and the number of islands satu- 
rates. 

Finally, we histogram the island sizes from all the in- 
dependent runs. Fig. [5] plots the average number of 
islands on the substrate against island size. For both 
models, a peak island size emerges for h > 2.5 MLs. 
For unfaceted islands, the histogram broadens signifi- 
cantly upon growth due to a wide distribution of growth 
rates. In contrast for faceted islands, it broadens much 
more slowly due to the highly kinetic limited growth 
mode. Nevertheless, the faceted islands do not possess 
narrower size distribution relative to the average size. 
This is because a significant size distribution already ex- 
ists when the islands become faceted as can be observed 
in Fig. m^b). The continuous nucleation of new islands 
also broadens the distribution as the older islands are 
larger on average. Another difference between the mod- 
els is that the peak of the histogram decays monotonically 
upon deposition for unfaceted islands while it increases 
for 2.5 < h < 4.5 due to the continuous nucleation of 
islands. 



V. DISCUSSIONS 

We have generalized a lattice model for strained films 
to allow for a range of local deformation states of sur- 
face atoms representing effective surface reconstructions. 
This deformations are assumed to be independent of the 
misfit induced strains for simplicity. Using this multi- 
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state lattice model, we have performed kinetic Monte 
Carlo simulations in 2D and observed the formation of 
(105) faceted pyramid islands. The model enables us to 
simulate faceted island formation in the kinetic limited 
regime. In this regime, island growth slows down dra- 
matically and becomes intermittent after faceting. The 
slower growth of the more established islands also leads 
to a continuous nucleation of islands until the substrate 
is fully occupied. The width of the island distribution 
is dominated both by fluctuations in the initial size at 
the start of faceting as well as the diversity in their ages. 
Stepped mounds from the conventional model exhibit a 
simple nucleation period followed by slow ripening. 

Additional studies on the growth and annealing of 
faceted islands under other growth conditions will be re- 
port elsewhere. It is also interesting to further generalize 
the model to consider two facet types so as to study the 
pyramid to dome transition. Generalization to 3D is con- 
ceptually straightforward but is challenging in practice 
because of the heavy computational load expected. 

This work was supported by HK RGC, Grant No. 
PolyU-5009/06P. 




FIG. 5: Size histograms for unfaceted (a) and faceted (b) 
islands 



[1] P. Politi, G. Grenet, A. Marty, A. Ponchet, and J. Villain, 
Phys. Rep. 324, 271 (2000). 

[2] V. A. Shchukin, N. N. Ledentsov, and D. Bimberg, Epi- 
taxy of nanostructures (Springer, 2003). 

[3] I. Berbezier and A. Ronda, Surf. Sci. Rep. 64, 47 (2009). 

[4] Y.-W. Mo, D. E. Savage, B. S. Swartzentruber, and M. G. 
Lagally, Phys. Rev. Lett. 65, 1020 (1990). 

[5] A. Vailionis, B. Cho, G. Glass, P. Desjardins, D. G. 
CahiU, and J. E. Greene, Phys. Rev. Lett. 85, 3672 
(2000). 

[6] M. Kastner and B. Voigtlander, Phys. Rev. Lett. 82, 2745 
(1999). 

[7] M. R. Mckay, J. A. Venables, and J. Drucker, Phys. Rev. 

Lett. 101 (2008). 
[8] A. Rastelh and H. von Kanel, Surf. Sci. 532, 769 (2003). 
[9] G. Medeiros-Ribeiro, A. Bratkovski, T. Kamins, 
D. Ohlberg, and R. Williams, Science 279, 353 (1998). 
[10] A. Rastelh, M. Stoffel, J. Tersofl, G. S. Kar, and O. G. 

Schmidt, Phys. Rev. Lett. 95, 026103 (2005). 
[11] Y. Fujikawa, K. Akiyama, T. Nagao, T. Sakurai, M. G. 
Lagally, T. Hashimoto, Y. Morikawa, and K. Terakura, 
Phys. Rev. Lett. 88, 176101 (2002). 
[12] S. Cereda, F. Montalenti, and L. Miglio, Surf. Sci. 591, 



23 (2005). 

[13] G.-H. Lu, M. Cuma, and F. Liu, Phys. Rev. B 72, 125415 
(2005). 

[14] O. Shklyaev, M. Beck, M. Asta, M. Miksis, and 

P. Voorhees, Phys. Rev. Lett. 94 (2005). 
[15] F. Montalenti, P. Raiteri, D. Migas, H. von Kanel, 

A. Rastelli, C. Manzano, G. Costantini, U. Denker, 

O. Schmidt, K. Kern, et al., Phys. Rev. Lett. 93 (2004). 
[16] S. Cereda and F. Montalenti, Phys. Rev. B 75, 195321 

(2007). 

[17] C. M. Retford, M. Asta, M. J. Miksis, P. W. Voorhees, 
and E. B. Webb III, Phys. Rev. B 75, 075311 (2007). 

[18] B. G. Orr, D. Kessler, C. W. Snyder, and L. Sander, 
Europhys. Lett. 19, 33 (1992). 

[19] K. E. Khor and S. Das Sarma, Phys. Rev. B 62, 16657 
(2000). 

[20] M. Meixner, E. SchoU, V. Shchukin, and D. Bimberg, 

Phys. Rev. Lett. 88 (2002). 
[21] C.-H. Lam, C.-K. Lee, and L. M. Sander, Phys. Rev. 

Lett. 89, 216102 (2002). 
[22] J. L. Gray, R. Hull, C.-H. Lam, P. Sutter, J. Means, and 

J. A. Floro, Phys. Rev. B 72, 155323 (2005). 
[23] M. T. Lung, C.-H. Lam, and L. M. Sander, Phys. Rev. 



7 



Lett. 95, 086102 (2005). 
[24] G. Russo and P. Smcrcka, .J. Comput. Phys. 214, 809 
(2006). 

[25] C.-H. Lam, M. T. Lung, and L. M. Sander, J. Sci. Com- 
put. 37, 73 (2008). 

[26] J. Y. Lee, M. J. Noordhoek, P. Smereka, H. McKay, and 
J. M. Millunchick, Nanotechnology 20 (2009). 

[27] R. Zhu, E. Pan, and P. W. Chung, Phys. Rev. B 75 



(2007). 

[28] C.-H. Lam and M. T. Lmig, Int. J. Mod. Phys. B 21, 
4219 (2007). 

[29] P. Liu, Y. Zhang, and C. Lu, Phys. Rev. B 68 (2003), 

ISSN 1098-0121. 
[30] J.-N. Aqua, T. Prisch, and A. Verga, Phys. Rev. B 76 

(2007). 



